
%% Constant rate process (Poisson Process)

T=100000;   %Number of small intervals
k=.01;      %Rate constant (s^-1)
dT=1;       %Time interval (s)
P=dT*k;     %Probability of event in each interval (P<<1)

%Time intevals are independent of each other, so we can do them all at once
Events=double(rand(T,1)<P)

%Visualize our events
figure
bar(Events)
xlabel('Time')
ylabel('Event')

%% Waiting time distribution

%find each event. Idx is the element index of non-zero values of Events
[Idx,Val]=find(Events)

%Find times in between each event
WaitTime=Idx(2:end)-Idx(1:end-1)

%Show histogram of wait times.
figure
hist(WaitTime)
xlabel('Time')
ylabel('Counts')

%% Counts distribution (Poisson Distribution)

%We have a different number of events each time we run the simulation
%Let's look at the distrubution of events. 

T=1000;         %Number of small time intervals
k=.05;          %Rate contant (s^-1)
dT=1;           %Time interval (s)
P=dT*k;         %Probability of event in each interval
lambda=T*P;     %Expected (average number) of events
Repeats=10000;  %Number of repeat trials of event sequence. 

Events=double(rand(T,Repeats)<P);

%find number of events for each repeat trial
Counts=sum(Events);

%Show the distribution
figure
hist(Counts,(0:100))
xlabel('Time')
ylabel('Counts')
s=sprintf('Lambda = %g',lambda)
title(s)

%% Generating Poisson Distributed Numbers

%Expected value
lambda=3;

%Vector with values equal to expected value
Vec=lambda*ones(10000,1);

%A vector of numbers selected from the Poisson distribution. 
PoisVec=poissrnd(Vec);

%Show Histogram
figure
hist(PoisVec,(0:100))
xlabel('Time')
ylabel('Counts')
s=sprintf('Lambda = %g',lambda)
title(s)

%% Generating Normally Distributed Numbers

mu=3;       %Expected value
sigma=10;   %Standard Deviation

%Vector with values equal to expected value
Vec=mu*ones(10000,1);

%A Vector of Nomrally distributed numbers
NormVec=normrnd(Vec,sigma);

%This is the same as:
%NormVec=Vec+sigma*randn(10000,1);

%Show Distribution
figure
hist(NormVec,100)
xlabel('Time')
ylabel('Counts')
s=sprintf('mu = %g; sigma = %g',mu,sigma)
title(s)

%% Generating Exponentially Distributed Numbers

%Expected value
lambda=3;

%Vector with values equal to expected value
Vec=lambda*ones(10000,1);

%Vector of Exponentially distributed values
ExpVec=exprnd(Vec);

%Look at distribution
figure
hist(ExpVec,100)
xlabel('Time')
ylabel('Counts')
s=sprintf('Lambda = %g',lambda)
title(s)

%% Calculating the values of the Probability Distribution Function (PDF)

%Poisson
lambda=3;   %Expected value
k=4;        %Test value
poisspdf(k,lambda)

%Exponential
lambda=3;   %Expected value
x=4;        %Test value
exppdf(x,lambda)

%Normal
mu=3;       %Expected value
sigma=10;    %Standard deviation of distribution
x=5;        %Test value
normpdf(x,mu,sigma)


%% Making Noisy Data from each distribution

%1st order polynomial
a=1;
b=5;

x=(0:100);
y=a+b*x;

%Poisson noise
y_P=poissrnd(y);

%Gaussian (normal) noise
sigma=5;
y_N=normrnd(y,sigma);

%Exponential
y_E=exprnd(y);

%Plot them all using 'subplot'
figure
subplot(2,2,1)
plot(x,y)
xlabel('x')
ylabel('y')
title('Expected')

subplot(2,2,2)
plot(x,y_P,'.')
xlabel('x')
ylabel('y')
title('Poisson')

subplot(2,2,3)
plot(x,y_N,'.')
xlabel('x')
ylabel('y')
title('Normal')

subplot(2,2,4)
plot(x,y_E,'.')
xlabel('x')
ylabel('y')
title('Exponential')





















